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Abstract 

This paper represents an overview of the most promising numerical methods for mul¬ 
tiphase modelling and its comparison on the example of known multiphase problems. 
Namely we will discuss Baer-Nunziato(BN) and recently proposed Godunov-Romosky- 
Peshkov(GPR) approaches. A distinctive and very important property of these ap¬ 
proaches is its hyperbolicity that provides a well-posedeness of the problem. In turn, 
it allows to apply well-developed numerical schemes. We will see that this property 
is achieved by the methods described by different ways, employing different ideas of 
formulation of governing equations. A complementary conservation equation is used in 
BN model, while GPR model takes an idea of applying of a distortion tensor to build 
a governing system of equations. It makes possible to obtain a unified mathematical 
model both for fluids aTd solids problem in frame of GPR model. That fact opens a 
very wide range of its application. 
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Introduction 


Multiphase flows are the omnipresent phenomena in our life. From particle interactions 
in molecular clusters up to behavior of star nebulas, in a number of scientific and engineer¬ 
ing disciplines. Ground processes in geology, airflows in aerospace tasks, reactive media in 
chemistry and biology, different scales of problems in astronomy - are just several examples 
of multiphase presence. A development of computational simulation of these problems, as 
an efficient tool to study different sides of the problem where a physical experiment can be 
intricate or not possible, is on behalf of many analytical and industrial research. 

Among all the methods proposed for the multiphase flows modelling one can pick out 
two remarkable ones that describe the flow via the hyperbolic system of equations. These are 
proposed by Baer and Nunziato model ([1], further denonted as BN) and model by Godunov 
and Romenski ([2]), recently advanced to the unified hyperbolic formulation of fluid and solid 
dynamics, in [3, 4] (further denoted as GPR model). 

Hyperbolicity of the system supplies the well-posedness of the problem, that is extremely 
important for mathematical modelling ([5]). 

In frame of BN model additional conservation relations are introduced to close the system 
and it gave quite accurate results for mixture, allowed to capture an interface in multiphase 
problems. BN model gives a hyperbolic system but not all of the equations are in conserva¬ 
tive form, it put some difficulties in the definition of the discontinuous solution and in the 
development of very high order numerical methods. 

GPR model presents both features- hyperbolicity and a conservative form of all the equa¬ 
tions of the system. These two specialties provide very good mathematical properties in¬ 
cluding continuous dependence of the solution on the initial data and the finite velocity for 
perturbation propagation. That makes model advantageous as it allows to apply known ac¬ 
curate numerical methods (like recent ADER scheme for example, [6, 7]). 

This paper represents a comparison between well-established BN and new GPR model 
approaches. In the sections below there are a brief description of both methods and com¬ 
parison of numerical experiments performed after these methods, a conclusion is given in the 
last part. 


Brief description of hyperbolic models 

Baer-Nunziato approach 

Several approaches were developed for numerical simulation of systems containing several 
different phases. Some systems of governing equations were constructed using an averaging 
of flow variables [8], using continuum mixture theory [9] and etc. In general all they revealed 
a closure problem for system of equation, as a supplementary independent kinematic variable 
(as volume fraction for example) is included. 
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To bypass that obstacle some additional conditions were applied, like relation for pressure 
equilibrium between the phases, or a definition pressure in the solid phase in terms of the 
configuration changes of the granular material (particularly for detonation problems). That 
imposed some constraints to the range of applicability of the methods. Substantially the 
difficulty of modelling was stuck on the issue of describing a compressibility of each phase 
and compaction of the solid phase. 

BN approach accounts for these effects using an additional conservation equation for the 
volume fraction of the solid with the help of evolutionary equation for the volume fraction 
consistent with thermodynamics, ft was firstly shown for the numerical simulation of the two 
phase reactive flows where a perturbation passed through the mixture of gas and small solid 
particles in frame of the deflagration or detonation processesQl]). The governing equations 
are the composition of two separate systems of the conservation laws for each phase, coupled 
by interface exchange terms ([10]). 

Thus obtained system of equations became hyperbolic. But not all the equations have a 
conservative form. For two phase flow system reads: 
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Here index m — 1,2- corresponds to phase 1 and 2 (ml and m2 - when there is a differ¬ 
ence between parameters of each phase), 5 m is equal to ’1’ for the first and ’-1’ for the second 
phase. a,p,u - are the volume fraction, density and velocity of each phase, p and A - are 
the relaxation parameters conformed to dynamic compaction viscosity and drag force (for 
pressure and velocity relaxation respectively). They are finite within two-phase flow regions 
and while instanteneous relaxation, and infinite at the interface. T), V); - pressure and velocity 
at the interface. 


GPR approach 


GPR model takes an idea of continuum particles (or material elements) of small finite 
scale as a means to describe a changing of the media. These continuum particles represent 
deformable ensembles of molecules existing during a finite time and rearrangements of which 
actually defines a flow phenomenon, [4], Elastic stretching of microscopic bonds between 


3 
















continuum particles results in internal resistance to flow (i.e. viscosity). The time of rear¬ 
rangement embodies a fluidity of a media. An admissible length scale and rearrangement 
time were proposed in [4, 11], Roughly saying, a large time values maps a solid material, 
time value tending to zero maps an ideal gas. 

Mathematical analysis [5| leads to a first order hyperbolic system of equations where vis¬ 
cous stresses are computed from the so-called distortion tensor A. That tensor A reflects 
a deformation of particles and is one of the primary state variables. A stiff source term in 
governing system accounts for strain relaxation in the evolution equations of A and rear¬ 
rangement time r. 

Rearrangement time r can be taken as constant for the Newtonian flows while the depen¬ 
dence of r on tensor A reflects the non-Newtonian properties of fluids transition from elastic 
to plastic regime in solids. 

The system of governing equations reads: 
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Where p, u, p, E, a retain their usual meaning of density, velocity, pressure, energy and 
stress tensor respectively. A - distortion tensor, source term i\)k m = ■ 9(r) - positive 

scalar function corresponding to the properties of studied material, that will be specified 
below. 


An important role in the system plays an energy definition, ft is suggested to be com¬ 
posed of three scale components. Microscale of medium’s molecules that depends on the 
local density and entropy (internal energy); mesoscale corresponding to the scale of material 
element (ensemble of molecules) - this one depends on the distortion tensor, and macroscale 
component for the flow scale (conventional kinetic energy). 


Equation of state employed in this work are as follows: 
- Internal energy (microscale part) Ei(p,s) 


Ei(p,s) 


7(7 - 1) 
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where q = 0 in case of ideal gas, and q = 1 in case of stiffened gas or liquid, 
indexes ’0’ correspond to some reference values 
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Mesoscale part E 2 (A) 


E 2 (s) = ^dev(G)dev(G) 
with dev(G) - trace-free part of the tensor dev( G) = A 7 A. 


The c s is the characteristic velocity of propagation of transverse perturbations. It may 
be laid as a constant approximate value or calculated from the value of phase velocity V (cj) 
corresponding to disappearance of dispersion, V(u) = \JCq + 4/3 c 2 s (an effect called as well 
as ’high frequency limit’). More detailed description can be found in [4, 3]. 

- Macroscale part, specific kinetic energy per unit mass E :i (u) 

E%[vj) — UiUi 

The definition of c s noted above is important as well for estimation of relaxation time r, 
standing in the right part of PDEs. The fact of the mater that relaxation time, or time of 
rearrangement of material elements, is a mark of the fluidity behaviour of the medium, hence 
a relation between time r and the viscosity coefficient rj takes place: r) = (see details 

in [3], section 3.1). 

Thus, deformable material elements with the ‘life time’ (time before rearrangement) defin¬ 
ing a fluidity of the medium permit to introduce a distortion tensor, that in its turn, gives 
a way to treat viscous fluids as generalized visco-plastic solids. This is a key aspect of GPR 
model. 


Numerical experiments 

We show in this section the results of GPR numerical solution of the one-dimensional 
shock-tube test with two fluids and its comparison with the results obtained in frame of BN 
model ([12]). 

The computational domain represents a shock tube filled on the left side with high- 
pressure liquid water and with air at normal conditions on the right side. Initial data are as 
follows: 

Left side: (water) p=1000kg/m 3 , p = 10 9 Pa, u =0 m/s, 7=4.4. 

Right side: (air) p=50kg/m 3 , p = 10 5 Pa, u =0 m/s, 7=1.4. 

A regular mesh with 1002 cells is used. 

The time step size is calculated via CFL number, which is kept with a value le-2 at the first 
10 steps of the calculation and than 0.9 for the rest. 

The master system (2) should be complemented with closure relations. For that purpose 
an equation of state is defined in the form of the generalized internal energy E as a function 
of parameters of state, described in a section above. Then pressure of each phase is computed 
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as a function of total energy derivative ad is as follows: : 



BE 

dp 


A high-order MUSCL-Hancock finite volume scheme in conjunction with path-conservative 
is employed. 


Graphics of phases pressure, phases velocity, volume fraction and mixture density for the 
time moment 229 ps are shown in the figure 1. The curves in red correspond to a GPR 
approach, blue lines correspond to BN model. A black dotted line displays an exact solution. 


The graphs of pressures, velocities and densities show in fact five curves: red curve, cor¬ 
responding to GPR model actually consists of two superposed curves (one for liquid and one 
for gas component), and blue curve as well consists of two curves of different components. 
The fifth curve (in black) is a single exact solution. In both GPR and BN cases these numer¬ 
ical curves for different phases are indiscernible. And that tell us a distinctive result of the 
pressure and velocity relaxation procedures, which give an accurate relaxed state. The phase 
volume fraction obtained via GPR model is compared with BN model and is not compared 
with the exact solution because no exact solution exists for the variables. 


One can see a good agreement between the results for phase pressures and velocities and 
their congruence with an exact solution. On the graphic of gas volume fraction we can note 
a slight difference of the perturbation propagation speed for GPR and BN model. A bit 
higher value of that speed for BN model than for GPR. That is due to the fact that finite 
wave speeds are included in GPR model by the definition of the last and their value are 
not affected by the numerical scheme effects, that in general cases is calculated taking into 
account different current variables meaning. 

An interesting scheme effect can be seen at the mixture density plot. A red (GPR) curve 
makes quite strong jump at the contact break up zone. This may make one conclude that 
a bug in a scheme takes place. However it becomes clear if we have a look at the density 
definition in GPR model,in system 2. The total density is calculated as a sum of weighted 
densities of its components ( ptotai = J2 a iPii where index i refers to a distinct component). 
Due to the fact that at the moment of rupture both volume concentrations are nonzero 
(according to the model definition with the distortion tensor presenting in all the domain), 
the density of the system takes on a value higher than elsewhere in the flow, therefore, a 
jump occurs. Nonetheless that effect has no influence on the other parameters of the flow. 
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Figure 1: Comparison of BN (in blue) and GPR (in red) model. Phase pressure (upper 
left-side picture), phase velocity (upper right-side picture), mixture density (bottom left-side 
picture), volume fraction (bottom right-side picture) 


Conclusion 

A brief but quite complete description of the most budding methods for multiphase mod¬ 
elling, namely BN and GPR, are presented in this paper. The numerical test of the same 
problem for both methods is demonstrated. Both methods describe the flows via hyperbolic 
system of equation. The hyperbolicity is a very important property that provides the well- 
posedeness of the problem. That, in its turn, leads to the opportunity to apply well-developed 
numerical schemes, their stability, accuracy and other useful features. 

BN model employs an additional conservation equation to solve a system closure problem 
that is regular for general multiphase models. It allows to find a solution for non-conservative 
terms, deals with the cases of strong shock waves or the flows with non equilibrium effects. 
This model applies to pure fluids and to mixtures. GPR model takes an idea of applying 
of a distortion tensor to build a system of governing equations. That makes possible to 
represent a unified mathematical model for both fluids and solids mechanics problems. An¬ 
other distinctive point of that model is its suitability for the use of advanced high-accuracy 
Godunov-type numerical methods. It gives an accurate solution for the main characteristics 
of multiphase problem. 

With the help of the demonstrated numerical test we can observe a good agreement of the 
GPR results with the ones of BN model, that reassures in its applicability and efficiency. For 
some non-conservative terms, as a volume fraction for example, GPR gives even more correct 
distribution and shear sound speed. Taking into account an accuracy of GPR model and 
its great area of application not just for fluids but for both liquid and solid media, authors 
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consider a new GPR approach a very promising approach and efficient tool for industrial and 
scientific tasks. 
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